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A  compact  empirical  mode  decomposition  (CEMD)  is  presented  to  reduce  mode  mixing, 
end  effect,  and  detrend  uncertainty  in  analysis  of  time  series  (with  N  data  points).  This 
new  approach  consists  of  two  parts:  (a)  highest-frequency  sampling  (HFS)  to  generate 
pseudo  extrema  for  effective  identification  of  upper  and  lower  envelopes,  and  (b)  a  set 
of  27V  algebraic  equations  for  determining  the  maximum  (minimum)  envelope  at  each 
decomposition  step.  Among  the  27V  algebraic  equations,  2(7V  —  2)  equations  are  derived 
on  the  base  of  the  compact  difference  concepts  using  the  Hermitan  polynomials  with  the 
values  and  first  derivatives  at  the  (TV  —  2)  non-end  points.  At  each  end  point,  zero  third 
derivative  and  determination  of  the  first  derivative  from  several  (odd  number)  nearest 
original  and  pseudo  extrema  provide  two  extra  algebraic  equations  for  the  value  and 
first  derivative  at  that  end  point.  With  this  well-posed  mathematical  system,  one  can 
reduce  the  mode  mixing,  end  effect,  and  detrend  uncertainty  drastically,  and  separate 
scales  naturally  without  any  a  priori  subjective  criterion  selection. 

Keywords :  Compact  empirical  mode  decompositions  (CEMD);  empirical  mode  decom¬ 
position  (EMD);  highest-frequency  sampling  (HFS);  pseudo  extrema;  compact  differ¬ 
ence;  highest-frequency  sampling  (HFS);  Hermitian  polynomials;  intrinsic  mode  function 
(IMF);  end  effect;  detrend  uncertainty;  mode  mixing. 


1.  Introduction 

Most  of  the  time  series  we  encounter  are  nonlinear  and  nonstationary.  Because  of 
the  limitations  of  available  methodologies  for  analyzing  data,  the  crucial  phase  of 
data  analysis  has  in  the  past  been  relegated  to  “data  processing”,  where  data  are 
routinely  put  through  certain  well-established  algorithms  to  extract  some  standard 
parameters.  These  traditional  methods  usually  assume  linearity  and  stationarity  of 
the  time  series.  Since  the  seminal  paper  by  Huang  et  al.  (1998),  adaptive  empirical 
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model  decomposition  (EMD)  has  been  developed  to  analyzing  nonlinear  and  non¬ 
stationary  data.  The  EMD  decomposes  a  nonlinear  and  nonstationary  signal  into 
several  intrinsic  mode  functions  (IMFs)  and  a  trend.  Instantaneous  frequency  can 
be  obtained  by  the  Hilbert-Huang  transform  (HHT),  and  then  the  time-frequency- 
energy  distribution  characteristics.  An  IMF  is  a  function  that  must  satisfy  two  con¬ 
ditions  according  to  the  algorithm  originally  developed:  (a)  the  difference  between 
the  number  of  local  extrema  and  the  number  of  zero-crossings  must  be  zero  or  one; 
(b)  the  running  mean  value  of  the  envelope  defined  by  the  local  maxima  and  the 
envelope  defined  by  the  local  minima  is  zero.  For  such  a  signal,  the  interior  extrema 
are  easily  identified.  However,  these  extrema  are  not  enough  to  determine  two  well- 
behaved  fitting  spline  envelopes  near  the  two  ends  for  the  sifting,  especially  in  the 
cases  when  the  total  number  of  splines  are  small,  for  the  extrapolation  of  a  spline 
often  leads  to  undesirable  big  error  especially  near  the  end  points. 

The  easiest  way  is  to  treat  the  two  end  points  as  “frozen”  points,  i.e.  the  two 
end  points  are  on  both  maximum  and  minimum  envelopes.  Such  a  treatment  makes 
the  trend  varying  from  the  first  end  point  to  the  last  end  point.  Other  practice  is 
to  extend  data  points  beyond  the  end  points  so  as  to  carry  out  the  spline  envelope 
fitting  over  and  even  beyond  the  existing  data  range;  such  as  the  wave  extension 
method  [Huang  et  al.  (1998)],  local  straight-line  extension  method  [Wu  and  Huang 
(2009)],  mirror  or  anti-mirror  extension  [Zhao  and  Huang  (2001)],  and  self-similarity 
[Wang  et  al.  (2007)],  to  name  a  few,  have  been  used.  While  methods  for  extending 
data  vary,  the  essence  of  all  these  methods  is  to  predict  data,  a  dauntingly  difficult 
procedure  even  for  linear  and  stationary  processes.  The  problem  that  must  be  faced 
is  how  to  make  predictions  for  nonlinear  and  nonstationary  stochastic  processes. 
The  original  signal  only  has  the  extrema  in  the  data  series,  so  the  extending  points 
beyond  the  two  end  points  are  not  real.  Therefore  the  data  extension  methods  will 
not  solve  the  problem  no  matter  how  much  efforts  have  been  spent.  Besides,  the 
end  error  may  propagate  from  the  ends  to  the  interior  of  the  data  span  that  would 
cause  severe  deterioration  of  the  IMFs  obtained. 

An  ensemble  EMD  was  proposed  to  solve  the  irregular  distribution  of  local 
maxima  (minima)  problem  [Wu  and  Huang  (2009)] .  This  approach  consists  of  sifting 
an  ensemble  of  white  noise-added  signal  (data)  and  treats  the  mean  as  the  final  true 
result.  The  reason  to  add  finite,  not  infinitesimal,  amplitude  white  noise  is  to  obtain 
relatively  uniform  distribution  of  local  maxima  and  minima  such  that  the  upper 
and  lower  envelope  can  be  identified. 

Questions  arise'.  Can  relatively  uniform  distribution  of  local  maxima  and  minima 
be  established  at  each  EMD  step  without  adding  any  white  noise?  Can  the  end 
effect  be  eliminated  without  using  either  extrapolation  or  interpolation  (with  extra 
point  beyond  the  end  point)?  Can  the  upper  and  lower  envelopes  be  determined 
in  a  systematic  way?  Can  the  trend  be  determined  objectively?  These  problems 
will  be  answered  in  this  study  through  identification  of  highest-frequency  pseudo 
extrema  and  use  of  compact  difference  concepts  [Chu  and  Fan  (1998),  (1999)]. 
With  highest-frequency  pseudo  extrema,  relatively  uniform  distributions  of  local 
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(original  and  pseudo)  maxima  and  minima  are  found.  Determination  of  either  upper 
or  lower  envelope  becomes  to  solve  a  set  of  2 N  algebraic  equations  of  the  values 
and  first  derivatives  for  each  envelope  at  all  N  data  points  (including  the  two  end 
points)  using  the  Hermitian  polynomials.  This  method,  called  the  compact  EMD 
(CEMD),  shows  evident  improvement  in  data  analysis.  The  rest  of  the  paper  is 
organized  as  follows.  Section  2  introduces  the  classical  EMD.  Section  3  depicts 
procedure  of  CEMD.  Section  4  demonstrates  the  reduction  of  the  end-point  effect 
by  CEMD.  Section  5  shows  the  enhancement  of  detrend  capability  by  CEMD. 
Section  6  presents  the  conclusions. 


2.  EMD 

Let  x(t)  represent  the  time  series  with  fluctuations  on  various  time  scales  (Fig.  1). 
The  EMD  method  is  depicted  as  follows.  First,  the  local  minima  and  maxima  of 
the  signal  x(t)  are  identified.  Second,  the  local  maxima  are  connected  together  by  a 
cubic  spline  interpolation  (other  interpolations  are  also  possible),  forming  an  upper 
envelope  emax(f).  The  same  is  done  for  local  minima,  providing  a  lower  envelope 
emin (t).  Third,  the  mean  of  the  two  envelopes  are  calculated 

rni(f)  =  [ema  x(f)  +  emin(t)]/2.  (1) 

Fourth,  the  mean  is  subtracted  from  the  signal,  providing  the  local  detail 

hi(t)  =  x{t)  -  TOi(t),  (2) 


Envelopes  and  the  Mean  :  data 


Time  :  second 

Fig.  1.  Procedure  of  traditional  EMD. 
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which  is  then  considered  to  check  if  it  satisfies  the  above  two  conditions  to  be  an 
IMF.  If  yes,  it  is  considered  as  the  first  IMF  and  denoted 


ci(i)  =  hi(t). 

It  is  subtracted  from  the  original  signal  and  the  first  residual, 

n(t)  =  x(t)  -  ci(f), 


(3) 


(4) 


is  taken  as  the  new  series  in  step  1.  If  hi(t)  is  not  an  IMF,  a  procedure  called 
“sifting  process”  is  applied  as  many  times  as  necessary  to  obtain  an  IMF.  In  the 
sifting  process,  hi(t)  is  considered  as  the  new  data,  and  the  same  procedure  applies. 
The  IMFs  are  orthogonal,  or  almost  orthogonal  functions  (mutually  uncorrelated). 
This  method  does  not  require  stationarity  and  linearity  of  the  data  and  is  especially 
suitable  for  nonstationary  and  nonlinear  time  series  analysis. 

By  construction,  the  number  of  extrema  decreases  when  going  from  one  residual 
to  the  next;  the  above  algorithm  ends  when  the  residual  has  only  one  extrema,  or 
is  constant,  and  in  this  case  no  more  IMF  can  be  extracted;  the  complete  decompo¬ 
sition  is  then  achieved  in  a  finite  number  of  steps.  The  signal  x(t)  is  finally  written 
as  the  sum  of  mode  time  series  Ci(t)  and  the  residual  rm(t ): 


m 


(5) 


There  is  no  any  oscillation  (i.e.  nonexistence  of  both  maximum  and  minimum 
envelopes)  in  the  residue  which  should  represent  the  trend.  Obviously,  suc¬ 

cessfulness  of  the  EMD  depends  on  accurate  determination  of  upper  and  lower 
envelopes.  Three  long  recognized  difficulties  (end-point  effect,  mode  mixing,  and 
detrend  uncertainty)  are  related  to  each  other  and  caused  by  uncertain  identifi¬ 
cation  of  upper  and  lower  envelopes  due  to  (a)  uncertain  maximum  and  minimum 
values  at  the  end  points  and  (b)  irregular  distribution  of  local  maximum  (minimum) 
points. 

3.  CEMD 

A  new  CEMD  method  is  developed  on  the  base  of  (a)  compact  difference  concept 
[e.g.  Chu  and  Fan  (1998),  (1999)]  using  Hermitian  polynomials  and  (b)  identifica¬ 
tion  of  uniformly  distributed  pseudo  maximum  and  minimum  points  using  highest- 
frequency  sampling.  CEMD  determines  upper  and  lower  envelopes  objectively  and 
accurately  at  all  N  data  points  including  the  two  end  points. 

3.1.  Basic  algebraic  equations 

Let  {xi,i  =  1,2,..., iV}  be  the  original  time  series  with  corresponding  time 
instances  {U,i  =  1,2,..., TV}  and  two  end  points  (x\,Xn)-  Let  (a^.max\j  =  1, 
2, . . . ,  J)  and  (a:[.min^ ,  k  =  1,2,...,  K)  be  the  local  maxima  and  minima  with  J  =  K 
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or  differing  at  most  by  one;  and  { m ,  i  =  1,2,...,  }  and  {U,  i  =  1,2,...,  }  be 

the  corresponding  upper  and  lower  envelopes.  Here,  the  upper  envelope  is  taken 
as  an  example  for  illustration.  On  the  upper  envelope,  the  values  are  given  at  the 
local  maxima  and  unknown  at  the  other  time  steps  (Fig.  2).  For  simplicity 

without  loss  of  generality,  a  uniform  time  step  At  =  ti+ 1  —  t*  (z  =  1,2,...,  JV  —  1) 
is  assumed. 

Let  be  the  value  and  first  derivative  of  the  upper  envelope  at  the 

time  ti .  For  an  interval  tj+i],  (z  =  2,  3, . . . ,  N  —  1),  a  quartic  Hermitian 
polynomial, 

Pa(0  =  a0  +  ai£  +  a2£2  +  a3£3  +  a4£4,  £  =  t~U,  (6) 

is  used  to  determine  (zz,,zzQ  simultaneously  (Fig.  3)  with 

p4(-At)  =  m- 1,  p4(  0)  =  Ui ,  p4(At)  =  ui+ 1, 

(7) 

pi(-M)  =  ui_1,  p'4(At)  =  ui+1. 

Substitution  of  (6)  into  (7)  leads  to  the  expressions  for  the  coefficients  of  the  Her¬ 
mitian  polynomials  (6), 


QjQ  —  'Uii 


°1  =  2 


Ui+l  —  Mi- 1 


a2  = 


2A  t 

Ui- 1  -  2 m  +  Uj+i  _  1 

At2  _  2 


-  -K+i  +<_i), 


M»+l  ~  Uj 

2At 
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Fig.  3.  Along  the  upper  envelope,  the  coefficients  of  the  quartic  Hermitian  polynomial  at  the 
point  i ,  is  calculated  from  ixj_ i,  Ui ,  Ui+ 1,  ^_1?  and 


«  3 


a4 


1 

4At2 

1 

2A  t2 


Uj+l  —  Ui- 1 

At 


H+l 


f  K+i  -  u'i~  i 
\  2At 


Ui- 1 


—  2rtj  +  Ui+i 

At2 


(8) 


First  derivative  of  (6)  at  time  ti  (i.e.  £  =  0)  gives 

p'M  =  a\  =  u'i. 


(9) 


Substitution  of  a\  from  (8)  into  (9)  leads  to 


Ui 


2  —  1 


4 


2+1 


=  6 


^z+1  1 

2At  ' 


(10) 


Since  ut  is  located  at  the  upper  envelope,  it  may  be  or  may  not  be  a  local  maximum. 
If  Ui  is  a  local  maximum,  it  has  a  given  value, 


Ui  =  x 


(max) 


(Ha) 


If  Ui  is  not  a  local  maximum,  an  additional  relation  is  needed  for  Ui.  This  relation 
can  be  obtained  by  reducing  the  order  of  the  Hermitian  polynomial  from  quartic  to 
cubic  (i.e.  the  fourth  derivative  equals  zero),  or  letting  04  =  0  in  (6),  which  leads  to 


Ui —  1  2 Ui  Ui- (_i  ^2+1  Ui—i 

At  ~  2 


(lib) 


Here,  Eq.  (8)  is  used.  The  values  and  first  derivatives  {ui,u[)  for  all  the  interior 
points  (*  =  2, 3, . . . ,  N  —  1)  are  obtained  by  solving  2  x  ( N  —  2)  algebraic  equations 
expressions  (10),  and  (11a)  or  (lib)  with  the  given  values  of  local  maxima  x jmax^ 
Two  more  algebraic  equations  are  needed  for  each  end  point  (i  =  1,  or  N).  The 
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first  equation  is  obtained  from  the  assumption  that  the  third  derivative  vanishes  at 
the  two  end  points, 


u"'  =  0,  u%  =  0. 


(12) 


Setting  the  quartic  Hermitian  polynomial  (6)  at  £2  (i.e.  i  =  2),  the  third  derivative 
at  ti  is  given  by 


d3P4(— At) 

dD 


=  60,3  +  24a,4  x  (—A t)  =  0. 


Substitution  of  (8)  with  i  =  2  into  (13)  leads  to 
9mi  —  16«2  +  7m3 


At 


=  3  m3  —  5m;. 


Similarly,  at  t n  we  have 

7um~2  —  16ujv-i  +  9mjv 
At 


=  5  %  —  3m; 


N  —  oaN_2- 


(13) 


(14) 


(15) 


The  second  equation  for  each  end  point  is  obtained  from  the  assumption  that  the 
first  derivatives  (G 1,  Gn )  can  be  calculated  at  the  two  end  points, 


u'i  —  Gi,  u'N  —  Gn. 


(16) 


3.2.  Determination  of  Gi,  Gn 

An  odd  number  (L)  of  extrema  points  ( x |°xt^ ,  l  =  1,  2, . . . ,  L)  close  to  each  end  point 
are  used  to  determine  Gi  and  Gn  with  a  weighted  least  squares  (WLS)  regression. 
Taking  G\  as  an  example,  the  WLS  regression  is  represented  by 

ijext  )=C  +  G1ti.  (17) 

In  the  WLS  regression,  cases  with  greater  weights  contribute  more  to  the  fit  of  the 
regression  line.  The  result  is  that  the  estimated  coefficients  are  usually  very  close  to 
what  they  would  be  in  ordinary  least  square  regression,  but  under  WLS  regression 
their  standard  errors  are  smaller.  Figure  4  shows  the  example  with  5  extrema  points 
for  determining  Gi.  The  weight  is  0.5  at  the  first  and  last  extrema  points  and  1.0 
otherwise.  Thus,  Gi  is  a  function  of  (xjext\  l  =  1,2, ... ,  L), 

g 4  =  iM4ext),4ext\  •  •  •  ,4ext))-  (18) 

Similarly,  Gn  is  a  function  of  (x[ext\l  =  M,  M  —  1, . . . ,  M  —  L  +  1), 

Gn  =  Fn  ,  x{M- i  >  ■  •  • .  X<M-K+ 1)  •  (19) 

Solutions  of  the  set  of  2 N  algebraic  equations  (10),  (11a),  (lib),  (14),  (15),  and 
(16)  with  (Gi,Gjv)  are  calculated  by  (18)  and  (19)  give  2 N  values  of  (m,,m()  for  all 
the  points  (i  =  1,  2, . . . ,  N). 
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3.3.  Pseudo  extrema  identification 

The  data  with  uneven  distribution  of  local  maximum  (minimum)  points  is  illus¬ 
trated  in  Fig.  5.  The  data  has  its  fundamental  part  as  a  low-frequency  sinusoidal 
wave  with  unit  amplitude.  At  the  two  middle  crests  and  first  three  trouphes  of  the 
low-frequency  wave,  high-frequency  intermittent  oscillations  with  amplitude  of  0.1 
are  riding  on  the  fundamental  [Fig.  5(a)].  Local  maximum  and  minimum  points 
are  distributed  unevenly.  The  high  and  low  frequency  components  are  shown  in 
Figs.  5(b)  and  5(c).  This  signal  has  no  trend  [Fig.  5(d)].  To  solve  such  a  prob¬ 
lem,  Wu  and  Huang  (2009)  proposed  an  ensemble  empirical  mode  decomposition 
(EEMD)  with  sifting  an  ensemble  of  white  noise-added  signal  (data)  and  treating 
the  mean  as  the  final  true  result.  This  method  creates  extra  local  maxima  and 
minima  in  such  a  way  that  a  uniform  reference  frame  in  the  time-frequency  space 
is  provided. 

Different  from  EEMD,  unevenly  distributed  local  maximum  and  minimum 
points  can  be  solved  by  identification  of  pseudo  maximum  and  minimum  points 
at  each  step  of  the  EMD  analysis.  Let  8  be  the  minimum  value  of  the  two 
smallest  half  time  period  of  two  neighboring  maximum  points  (5i)  and  minimum 
points  (<52), 


S  =  min(i5i,  S2), 


(20) 


which  represents  highest  frequency  fluctuation.  For  a  neighboring  local  maximum- 
minimum  pair  away  from  the  end  points  with  A  the  time  period  in  between  (Fig.  6), 
the  integer  p  is  calculated  by 


A 

J 


(21) 
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Fig.  5.  Time  series  containing  a  high  frequency  component  and  a  low  frequency  component 
without  trend:  (a)  signal,  (b)  high  frequency  component,  (c)  low  frequency  component,  and  (d) 
zero  trend. 
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Fig.  6.  Pseudo  maxima  and  minima  generated  by  highest- frequency  sampling  for  the  upper 
envelope:  (a)  local  maximum  away  from  the  end  point,  and  (b)  local  maximum  next  to  the  end 
point.  It  is  noted  that  similar  procedure  is  adopted  if  the  local  minimum  on  the  left  side  of  the 
pair. 


where  the  bracket  [  ]  shows  the  integer  part  of  the  real  number  inside  the  bracket. 
Pseudo  maximum-minimum  pairs  (p-pairs)  are  identified  with  the  time  period, 


6  = 


A 

2p  +  1  ’ 


(22) 


between  the  two  neighboring  pseudo  maximum  and  minimum.  For  the  local  maxi¬ 
mum  or  minimum  next  to  the  end  points  (local  maximum  point  as  an  example  in 
Fig.  6),  the  integer  (q)  is  calculated  by 


q  = 


A 

J 


-  1. 


Pseudo  maximum  and  minimum  points  are  identified  between  local  maximum  (or 
minimum)  and  the  end  point  with  the  time  period, 


6  = 


A 

9  +  1’ 


(23) 


between  the  two  neighboring  pseudo  maximum  and  minimum. 
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Fig.  7.  Identification  of  pseudo  maximum  and  minimum  points  by  highest  frequency  sampling: 
(a)  original  signal  [same  as  in  Fig.  6(a)],  (b)  distribution  of  pseudo  maxima  (red  circles)  and 
minima  (blue  circles),  (c)  upper  (red)  and  lower  (blue)  envelopes,  and  (d)  sifting. 


The  identified  pseudo  extrema  have  the  following  features:  (a)  the  frequency 
of  the  pseudo  maximum  (or  minimum)  points  is  nearly  the  same  and  not  higher 
than  the  highest  frequency  of  original  signal;  (b)  the  total  (i.e.  original  and  pseudo) 
maximum  and  minimum  points  must  be  alternatively  distributed  [Fig.  7(a)];  (3) 
the  pseudo  maximum  and  minimum  still  keep  the  original  values;  and  (4)  the  total 
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(original  and  pseudo)  extrema  will  be  distributed  uniformly  [Fig.  7(b)].  With  the 
uniform  distributed  total  extrema,  the  upper  and  lower  envelopes  can  be  objectively 
identified  [Fig.  7(c)].  After  sifting,  the  high-frequency  component  is  easily  obtained 
[Fig.  7(d)], 


3.4.  Elimination  of  mode  mixing 

The  CEMD  is  conducted  on  the  data  shown  in  Fig.  5(a).  Two  IMFs  and  a  residue 
(zero  line)  are  obtained  (Fig.  8).  Comparison  between  Figs.  8  and  5  shows  the 
capability  of  the  CEMD  to  analyze  the  original  time  series,  f(t)  =  fi{t)  +  /2(f), 
shown  in  Fig.  5.  During  the  comparison,  the  correlation  coefficient  (CC)  and  the 
relative  root  mean  square  error  (RRMSE), 


RRMSE  = 


(24) 


are  used  with  and  being  the  time  series  of  components  (trend  and  compo- 
nents/IMFs)  of  original  and  analyzed  data,  N  the  number  of  data  points,  and 
the  mean  of  { T f0'*}. 


CEMD  (CPU  time=0.328s) 


CorrCoef:  0.921. 
RRMSE:  0.444 

1 

1  1  1 

200  400  600  800 

Time 


Fig.  8.  CEMD  analysis  on  the  time  series  {xi}  given  in  Fig.  5(a)  with  IMF-1  resembling  the  high- 
frequency  component  /i(t)  in  Fig.  5(b),  IMF-2  resembling  the  low-frequency  component  /2M,  and 
no  trend  and  IMF- 2  resembling.  It  clearly  shows  the  capability  of  eliminating  mode  mixing,  end 
effect,  and  detrend  uncertainty. 
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The  original  two  components  (/i  and  f%)  in  Fig.  5  are  fully  recovered  by  the 
CEMD.  The  IFM-1  [i.e.  ci(f)]  resembles  fi(t)  with  CC  of  0.92  and  RRMSE  of  0.44. 
The  IFM-2  [i.e.  C2{t)]  resembles  /2(f)  in  Fig.  5  with  CC  of  1.0  and  RRMSE  of  0.03. 
Moreover,  the  CEMD  provides  the  one-to-one  correspondence  of  the  components 
between  original  analyzed  data:  C\(t)  to  02(f)  to  /2(f),  zero  trend  to  zero  trend. 


4.  Trend  Determination 
4.1.  Signal  with  a  trend 

A  time  series  consisting  of  a  quadratic  trend  and  four  harmonics, 

3 

x{ti)  =  f0(ti)  +  Y  f0(U)  =  A0tf, 

k=  1 

fk(ti)  =  Aksm(ukti  +  <pk),  fe  =  1,2,3  (25) 

Xi=x(ti),  ti  =  (i  —  l)At,  ti  =  0, 

fjv  =  0.9s,  At  =  0.0018  s,  TV  =  501, 

is  used  to  demonstrate  the  end  effect  in  EMD.  The  parameters  in  (25)  are  given  in 
Table  1.  Figure  9  shows  the  time  series  of  {xt\  with  the  two  end  points  with  the 
values  of 


x\  =  0.1096,  xn  =  0.4495. 


(26) 


Table  1.  Values  of  parameters  used  in  Eq.  (25). 
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Fig.  9.  Time  series  {xi}  represented  by  Eq.  (25). 
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Figure  10  shows  the  trend  and  harmonics  of  the  data  shown  in  Fig.  9:  foi  =  fo(U), 
fu  =  fi(ti),  hi  =  h(U),  and  f3i  =  f3(U ).  Obviously,  only  f0(U)  represents  the 
trend  of  the  signal  {xi}.  The  trend  of  {x^}  varies  from  /0i  to  /oat, 

/or  =  0,  foN  =  0.405.  (27) 

The  CEMD  is  conducted  on  the  time  series  {xi}  (shown  in  Fig.  9)  to  obtain  three 
IMFs  and  a  trend  (Fig.  11).  Obviously,  IMFs  well  correspond  to  the  harmonics,  Ci(t) 
versus  /*(£),  with  high  CCs  (0.9992  between  ci  and  /i,  0.9836  between  C2  and  hi 
0.9992  between  c3  and  /3,  0.9890  between  two  trends)  and  low  RRMSEs  (0.1834 
between  C\  and  /i,  0.082  between  c-i  and  hi  0-0441  between  c3  and  f3,  0.172  between 
two  trends). 

4.2.  Existence  of  noises 

Comparison  between  trends  from  the  original  data  [i.e.  the  quadratic  trend  /o(ti)]> 
the  CEMD  (Fig.  11)  clearly  shows  perfect  detrend  using  CEMD  and  false  detrend 
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Fig.  10.  Components  of  time  series  f.r, }:  a  quadratic  trend  and  four  harmonics  given  by  Eq.  (25). 
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Fig.  11.  CEMD  on  the  signal  shown  in  Fig.  9:  (a)  original  time  series,  (b)  IMF-1,  (c)  IMF-2, 
(d)  IMF-3,  and  (e)  trend.  Comparison  between  Figs.  11  and  10  shows  the  capability  to  reduce  the 
mode  mixing,  end  effect,  and  detrend  uncertainty. 


using  the  traditional  EMD  and  EEMD.  Since  the  observational  data  contain 
errors,  four  time  series  sm(ti)  (m  =  1,2,3)  are  constructed  each  by  a  signal 
[components  of  (25)]  and  random  noise  n{ti)  produced  by  a  pseudorandom  number 
generator  with  the  amplitude  of  0.1, 

si{ti)  =  f0(ti)  +  fi(U)  +  f2(ti)  +  n(ti), 

S2  {U)  =  fo(U)  +  fi{U)  +  f2(U )  +  f3(ti)  +  n(U ),  (28) 

■S3  {U)  =  fo{ti)  +  fliti)  +  f2{ti)  +  f3(ti)  +  +  n(ti). 

Here,  the  time  step  and  the  total  number  of  data  (N)  is  the  same  as  shown  in  (25). 
Figure  12  shows  the  three  time  series,  sm(U)- 
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Fig.  12.  Time  series  of  three  noisy  signals  (si,  s2 .  S3). 


Figure  13  also  shows  the  capability  of  CEMD  in  detrend  with  high  positive  CCs, 
CCi  =  0.984,  CC2  =  0.989,  CC3  =  0.996,  (29) 

and  low  RRMSEs, 

RRMSEi  =  0.287,  RRMSE2  =  0.359,  RRMSE3  =  0.502.  (30) 

The  ratio  of  CPU  with  EEMD-lk  versus  CEMD  is  310  (77.5  s  versus  0.25  s)  for 
analyzing  si(t),  270  (80s  versus  0.296s)  for  analyzing  s2(f),  and  220  (79.1s  versus 
0.359s)  for  analyzing  s3(t).  The  ratio  of  CPU  with  EEMD-5k  versus  CEMD  is  1548 
(387s  versus  0.25  s)  for  analyzing  si(t),  1365  (404  s  versus  0.296  s)  for  analyzing 
s2(t),  and  1103  (396  s  versus  0.359s)  for  analyzing  s3(f). 
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Fig.  13.  Detrend  of  three  noisy  signals  (si,  S2 ,  S3)  by  CEMD.  Here,  in  each  panel,  the  bold 
dashed-dotted  curve  is  the  trend  of  signal,  /o(t).  The  solid  curve  is  the  trend  identified  using 
CEMD. 

5.  Conclusions 

(1)  Mode  mixing,  end  effect,  and  related  detrend  uncertainty  in  the  traditional  EMD 
and  EEMD  can  be  eliminated  using  the  new  method  (i.e.  CEMD),  which  has  been 
established  on  the  base  of  (a)  compact  difference  scheme  concepts  [Chu  and  Fan 
(1998),  (1999)]  and  (b)  identification  of  uniformly  distributed  pseudo  maximum  and 
minimum  points  using  highest-frequency  sampling.  Determination  of  either  upper 
or  lower  envelope  becomes  a  mathematical  problem  to  solve  a  set  of  2N  algebraic 
equations  of  the  values  and  first  derivatives  for  each  envelope  at  all  N  data  points 
(including  the  two  end  points). 
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(2)  The  CEMD  uses  hybrid  Hermitain  polynomials  rather  than  the  cubic  spline. 
A  set  of  2N  algebraic  equations  (m,  it')  for  the  upper  envelope  (same  for  the  lower 
envelope)  are  obtained  from  more  realistic  conditions  such  as  zero  third  derivative 
and  use  of  weighted  regressed  first  derivative  at  the  end  points.  No  extending  points 
are  needed. 

(3)  Capability  of  the  CEMD  for  eliminating  mode  mixing,  end  effect,  and 
detrend  uncertainty  is  demonstrated  using  a  time  series  consisting  of  a  quadratic 
trend  and  four  harmonics.  The  CEMD  can  recover  the  original  signal  really  well 
with  low  RRMSEs.  Furthermore,  the  capability  of  the  CEMD  is  not  affected  by  the 
existence  of  random  noises. 


References 

Chu,  P.  C.  and  Fan,  C.  W.  (1998).  A  three-point  combined  compact  difference  scheme.  J. 
Comput.  Phys.,  140:  370-399. 

Chu,  P.  C.  and  Fan,  C.  W.  (1999).  A  three-point  non-uniform  combined  compact  difference 
scheme.  J.  Comput.  Phys.,  148:  663-674. 

Huang,  N.  E,  Shen,  Z,  Long,  S.  R,  Wu,  M.  C.,  Shih,  H.  H,  Zheng,  Q.,  Yen,  N.-C,  Tung, 
C.  C.  and  Liu,  H.  H.  (1998).  The  empirical  mode  decomposition  and  the  Hilbert 
spectrum  for  nonlinear  and  nonstationary  time  series  analysis.  Proc.  R.  Soc.  London, 
Ser.  A,  454:  903-993. 

Wang,  J.,  Peng,  Y.  and  Peng,  X.  (2007).  Similarity  searching  based  boundary  effect  pro¬ 
cessing  method  for  empirical  mode  decomposition.  Elect.  Lett.,  43(1):  58-59. 

Wu,  Z.  and  Huang,  N.  E.  (2009).  Ensemble  empirical  mode  decomposition:  A  noise-assisted 
data  analysis  method.  Adv.  Adapt.  Data  Analys.,  1:  1-41. 

Zhao,  J.  and  Huang,  D.  (2001).  Mirror  extending  and  circular  spline  function  for  empirical 
mode  decomposition  method.  J.  Zhejiang  Univ.  Sci.,  2:  247-252. 


1250017-18 


